home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / pascal / tpl60n19.zip / TESTPRGS.ZIP / MACHARIT.PAS < prev    next >
Pascal/Delphi Source File  |  1993-02-14  |  14KB  |  368 lines
  1. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  2.  
  3. UNIT MachArit;   { converted from Fortran original 05-01-92 Norbert Juffa }
  4.  
  5. INTERFACE
  6.  
  7. PROCEDURE MACHAR (VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  8.                       MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
  9.  
  10. PROCEDURE PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  11.                       MAXEXP: LONGINT; EPS,EPSNEG,XMIN,XMAX: REAL);
  12.  
  13.  
  14. IMPLEMENTATION
  15.  
  16. PROCEDURE MACHAR(VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  17.                      MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
  18.  
  19.  
  20. {-----------------------------------------------------------------------
  21. C  This Fortran 77 subroutine is intended to determine the parameters
  22. C   of the floating-point arithmetic system specified below.  The
  23. C   determination of the first three uses an extension of an algorithm
  24. C   due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
  25. C   but not all, of the improvements suggested by M. Gentleman and S.
  26. C   Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
  27. C   program was published in the book Software Manual for the
  28. C   Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
  29. C   Englewood Cliffs, NJ, 1980.  The present version is documented in
  30. C   W. J. Cody, "MACHAR: A subroutine to dynamically determine machine
  31. C   parameters," TOMS 14, December, 1988.
  32. C
  33. C  The program as given here must be modified before compiling.  If
  34. C   a single (double) precision version is desired, change all
  35. C   occurrences of CS (CD) in columns 1 and 2 to blanks.
  36. C
  37. C  Parameter values reported are as follows:
  38. C
  39. C       IBETA   - the radix for the floating-point representation
  40. C       IT      - the number of base IBETA digits in the floating-point
  41. C                 significand
  42. C       IRND    - 0 if floating-point addition chops
  43. C                 1 if floating-point addition rounds, but not in the
  44. C                   IEEE style
  45. C                 2 if floating-point addition rounds in the IEEE style
  46. C                 3 if floating-point addition chops, and there is
  47. C                   partial underflow
  48. C                 4 if floating-point addition rounds, but not in the
  49. C                   IEEE style, and there is partial underflow
  50. C                 5 if floating-point addition rounds in the IEEE style,
  51. C                   and there is partial underflow
  52. C       NGRD    - the number of guard digits for multiplication with
  53. C                 truncating arithmetic.  It is
  54. C                 0 if floating-point arithmetic rounds, or if it
  55. C                   truncates and only  IT  base  IBETA digits
  56. C                   participate in the post-normalization shift of the
  57. C                   floating-point significand in multiplication;
  58. C                 1 if floating-point arithmetic truncates and more
  59. C                   than  IT  base  IBETA  digits participate in the
  60. C                   post-normalization shift of the floating-point
  61. C                   significand in multiplication.
  62. C       MACHEP  - the largest negative integer such that
  63. C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
  64. C                 MACHEP is bounded below by  -(IT+3)
  65. C       NEGEPS  - the largest negative integer such that
  66. C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
  67. C                 NEGEPS is bounded below by  -(IT+3)
  68. C       IEXP    - the number of bits (decimal places if IBETA = 10)
  69. C                 reserved for the representation of the exponent
  70. C                 (including the bias or sign) of a floating-point
  71. C                 number
  72. C       MINEXP  - the largest in magnitude negative integer such that
  73. C                 FLOAT(IBETA)**MINEXP is positive and normalized
  74. C       MAXEXP  - the smallest positive power of  BETA  that overflows
  75. C       EPS     - the smallest positive floating-point number such
  76. C                 that  1.0+EPS .NE. 1.0. In particular, if either
  77. C                 IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  78. C                 Otherwise,  EPS = (FLOAT(IBETA)**MACHEP)/2
  79. C       EPSNEG  - A small positive floating-point number such that
  80. C                 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
  81. C                 or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  82. C                 Otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
  83. C                 NEGEPS is bounded below by -(IT+3), EPSNEG may not
  84. C                 be the smallest number that can alter 1.0 by
  85. C                 subtraction.
  86. C       XMIN    - the smallest non-vanishing normalized floating-point
  87. C                 power of the radix, i.e.,  XMIN = FLOAT(IBETA)**MINEXP
  88. C       XMAX    - the largest finite floating-point number.  In
  89. C                 particular  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  90. C                 Note - on some machines  XMAX  will be only the
  91. C                 second, or perhaps third, largest number, being
  92. C                 too small by 1 or 2 units in the last digit of
  93. C                 the significand.
  94. C
  95. C     Latest revision - December 4, 1987
  96. C
  97. C     Author - W. J. Cody
  98. C              Argonne National Laboratory
  99. C
  100. C-----------------------------------------------------------------------}
  101.  
  102. VAR   I, L, ITEMP, IZ,
  103.       J, K, MX, NXRES:     LONGINT;
  104.       A, B, BETA, BETAIN,
  105.       BETAH, ONE, T, TEMP,
  106.       TEMPA, TEMP1, TWO,
  107.       Y, Z, ZERO:          REAL;
  108.       CONV:                ARRAY [0..10] OF REAL;
  109.  
  110. LABEL 1, 2, 10, 21, 22, 30, 32, 40,
  111.       41, 42, 43, 44, 45, 46, 50, 52;
  112.  
  113. BEGIN
  114.  
  115. {-----------------------------------------------------------------------}
  116.  
  117.    FOR L := 1 TO 10 DO BEGIN
  118.       CONV [L] := L;
  119.    END;
  120.  
  121.    ONE  := CONV [1];
  122.    TWO  := ONE + ONE;
  123.    ZERO := ONE - ONE;
  124.  
  125. {-----------------------------------------------------------------------}
  126. {  Determine IBETA, BETA ala Malcolm.                                   }
  127. {-----------------------------------------------------------------------}
  128.  
  129.    A := ONE;
  130. 1: A := A + A;
  131.    TEMP  := A+ONE;
  132.    TEMP1 := TEMP-A;
  133.    IF TEMP1-ONE = ZERO THEN
  134.       GOTO 1;
  135.    B := ONE;
  136. 2: B := B + B;
  137.    TEMP := A+B;
  138.    ITEMP := TRUNC (TEMP-A);
  139.    IF ITEMP = 0 THEN
  140.       GOTO 2;
  141.    IBETA := ITEMP;
  142.    BETA  := CONV [IBETA];
  143.  
  144. {-----------------------------------------------------------------------}
  145. {  Determine IT, IRND.                                                  }
  146. {-----------------------------------------------------------------------}
  147.  
  148.    IT := 0;
  149.    B  := ONE;
  150. 10:IT := IT + 1;
  151.    B := B * BETA;
  152.    TEMP := B+ONE;
  153.    TEMP1 := TEMP-B;
  154.    IF TEMP1-ONE = ZERO THEN
  155.       GOTO 10;
  156.    IRND := 0;
  157.    BETAH := BETA / TWO;
  158.    TEMP := A+BETAH;
  159.    IF TEMP-A <> ZERO THEN
  160.       IRND := 1;
  161.    TEMPA := A + BETA;
  162.    TEMP := TEMPA+BETAH;
  163.    IF (IRND = 0) AND (TEMP-TEMPA <> ZERO) THEN
  164.       IRND := 2;
  165.  
  166. {-----------------------------------------------------------------------}
  167. {  Determine NEGEP, EPSNEG.                                             }
  168. {-----------------------------------------------------------------------}
  169.  
  170.    NEGEP := IT + 3;
  171.    BETAIN := ONE / BETA;
  172.    A := ONE;
  173.    FOR I := 1 TO NEGEP DO BEGIN
  174.       A := A * BETAIN;
  175.    END;
  176.    B := A;
  177. 21:TEMP := ONE-A;
  178.    IF TEMP-ONE <> ZERO THEN
  179.       GOTO 22;
  180.    A := A * BETA;
  181.    NEGEP := NEGEP - 1;
  182.    GOTO 21;
  183. 22:NEGEP := -NEGEP;
  184.    EPSNEG := A;
  185.  
  186. {-----------------------------------------------------------------------}
  187. {  Determine MACHEP, EPS.                                               }
  188. {-----------------------------------------------------------------------}
  189.  
  190.    MACHEP := -IT - 3;
  191.    A := B;
  192. 30:TEMP := ONE+A;
  193.    IF TEMP-ONE <> ZERO THEN
  194.       GOTO 32;
  195.    A := A * BETA;
  196.    MACHEP := MACHEP + 1;
  197.    GOTO 30;
  198. 32:EPS := A;
  199.  
  200. {-----------------------------------------------------------------------}
  201. {  Determine NGRD.                                                      }
  202. {-----------------------------------------------------------------------}
  203.  
  204.    NGRD := 0;
  205.    TEMP := ONE+EPS;
  206.    IF (IRND = 0) AND (TEMP*ONE-ONE <> ZERO) THEN
  207.       NGRD := 1;
  208.  
  209. {-----------------